home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MATH.SWG / 0033_CALCULUS.pas < prev    next >
Pascal/Delphi Source File  |  1993-11-02  |  5KB  |  155 lines
  1. { Updated NUMBERS.SWG on November 2, 1993 }
  2.  
  3. {
  4. LOU DUCHEZ
  5.  
  6. Hey everybody!  This unit performs calculus operations via basic numerical
  7. methods : integrals, derivatives, and extrema.  By Lou DuChez.  I don't
  8. want any money for this; please just leave my name in the source code
  9. somewhere, since this is the closest I'll ever get to being famous.
  10.  
  11. All functions return real values.  The last parameter in each function is
  12. a pointer to a "real" function that takes a single "real" parameter:
  13. for example, y(x).  See prior message to Timothy C. Novak for sample prog }
  14.  
  15. unit calculus;
  16. interface
  17.  
  18. function integral(a, b, h : real; f : pointer) : real;
  19. function derivative(x, dx : real; f : pointer) : real;
  20. function extremum(x, dx, tolerance : real; f : pointer) : real;
  21.  
  22. implementation
  23.  
  24. type
  25.   fofx = function(x : real) : real;     { needed for function-evaluating }
  26.  
  27. function integral(a, b, h : real; f : pointer) : real;
  28. var
  29.   x, summation : real;
  30.   y            : fofx;
  31. begin                                 { Integrates function from a to b,  }
  32.   @y := f;                            { by approximating function with    }
  33.   summation := 0;                     { rectangles of width h. }
  34.   x := a + h/2;
  35.   while x < b do
  36.   begin                               { Answer is sum of rectangle areas, }
  37.     summation := summation + h*y(x);  { each area being h*y(x).  X is at  }
  38.     x := x + h;                       { the middle of the rectangle.      }
  39.   end;
  40.   integral := summation;
  41. end;
  42.  
  43. function derivative(x, dx : real; f : pointer) : real;
  44. var
  45.   y : fofx;
  46. begin                 { Derivative of function at x: delta y over delta x }
  47.   @y := f;                                       { You supply x & delta x }
  48.   derivative := (y(x + dx/2) - y(x - dx/2)) / dx;
  49. end;
  50.  
  51.  
  52. function extremum(x, dx, tolerance : real; f : pointer) : real;
  53. { This function uses DuChez's Method for finding extrema of a function (yes,
  54.   I seem to have invented it): taking three points, finding the parabola
  55.   that connects them, and hoping that an extremum of the function is near
  56.   the vertex of the parabola.  If not, at least you have a new "x" to try...
  57.  
  58.   X is the initial value to go extremum-hunting at; dx is how far on either
  59.   side of x to look.  "Tolerance" is a parameter: if two consecutive
  60.   iterations provide x-values within "tolerance" of each other, the answer
  61.   is the average of the two. }
  62. var
  63.   y           : fofx;
  64.   gotanswer,
  65.   increasing,
  66.   decreasing  : boolean;
  67.   oldx        : real;
  68.   itercnt     : word;
  69. begin
  70.   @y := f;
  71.   gotanswer := false;
  72.   increasing := false;
  73.   decreasing := false;
  74.   itercnt := 1;
  75.   repeat                               { repeat until you have answer }
  76.     oldx := x;
  77.     x := oldx - dx*(y(x+dx) - y(x-dx)) /    { this monster is the new value }
  78.          (2*(y(x+dx) - 2*y(x) + y(x-dx)));  { of "x" based DuChez's Method }
  79.     if abs(x - oldx) <= tolerance then
  80.       gotanswer := true                     { within tolerance: got an answer }
  81.     else
  82.     if (x > oldx) then
  83.     begin
  84.       if decreasing then
  85.       begin              { If "x" is increasing but it }
  86.         decreasing := false;                { had been decreasing, we're }
  87.         dx := dx/2;                         { oscillating around the answer. }
  88.       end;                                { Cut "dx" in half to home in on }
  89.       increasing := true;                   { the extremum. }
  90.     end
  91.     else
  92.     if (x < oldx) then
  93.     begin
  94.       if increasing then
  95.       begin              { same thing here, except "x" }
  96.         increasing := false;                { is now decreasing but had }
  97.         dx := dx/2;                         { been increasing }
  98.       end;
  99.       decreasing := true;
  100.     end;
  101.   until gotanswer;
  102.  
  103.   extremum := (x + oldx) / 2;               { spit out answer }
  104. end;
  105.  
  106. end.
  107.  
  108.  
  109.  
  110. {
  111. I've put together a unit that does calculus.  This unit could be used, for
  112. example, to approximate the area under a curve (like a circle).
  113.  
  114. Because of the funny way my offline reader breaks up messages, I'm going
  115. to send you a "test" program first -- which just happens to calculate
  116. the area under a quarter circle -- then the following message (I hope)
  117. will be the unit source code.
  118. }
  119.  
  120. program mathtest;
  121. uses
  122.   calculus;
  123.  
  124. var
  125.   answer : real;
  126.  
  127. {$F+}                       { WARNING!  YOU NEED "FAR" FUNCTIONS! }
  128. function y(x : real) : real;
  129. begin
  130.   y := 4 * sqrt(1 - x * x);
  131. end;
  132.  
  133. begin
  134.   writeln('Function: y = (1 - x^2)^(1/2) (i.e., top half of a circle)');
  135.   writeln;
  136.  
  137. { Calc operations here are: }
  138.  
  139. { Integrate function from 0 to 1, in increments of 0.001. A quarter circle. }
  140. { Get slope of function at 0 by evaluating points 0.01 away from each other. }
  141. { Find extremum of function, starting at 0.4, initially looking at points
  142.   0.1 on either side of 0.4, and not stopping until we have two x-values
  143.   within 0.001 of each other. }
  144.  
  145.   answer := integral(0, 1, 0.001, @y);
  146.   writeln('Integ: ', answer:13:9);
  147.  
  148.   answer := derivative (0, 0.01, @y);
  149.   writeln('Deriv: ', answer:13:9);
  150.  
  151.   answer := extremum(0.4, 0.1, 0.001, @y);
  152.   writeln('Extrm: ', answer:13:9);
  153. end.
  154.  
  155.